topo_change <- st_read('data/topo_change/topo_change_polygons.shp')
## Reading layer `topo_change_polygons' from data source `C:\Users\keman\Dropbox\Documents\Courses\WR674\Mining Impacts\mining_impacts_to_elevation\data\topo_change\topo_change_polygons.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8958 features and 51 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -2332814 ymin: 407783 xmax: 2238422 ymax: 3128340
## epsg (SRID): NA
## proj4string: +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
#Checkout the column names
#names(topo_change)
biggest_area <- topo_change %>%
arrange(desc(AREA_SQ_KM)) %>%
slice(1:10)
# Map all mines
mapview(biggest_area)
# Subset to just Arizona
az_mine <- biggest_area %>%
filter(QUADNAME == 'esperanza_mill_AZ')
#Check that it is the right site
mapview(az_mine)
#Check projection of az_mine
#st_crs(az_mine)
az_raster_before <- get_elev_raster(az_mine,z=12)
## Merging DEMs
## Reprojecting DEM to original projection
## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
## +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
#Look at the structure of the data
#str(az_raster)
#Summary of the data
#summary(az_raster)
plot(az_raster_before)
#mapview(az_raster_before)
small_az <- aggregate(az_raster_before,fact=6)
mapview(small_az)
az_slope <- terrain(small_az,opt='slope',unit='degrees')
mapview(az_slope)
az_aspect <- terrain(small_az,opt='aspect',unit='degrees')
mapview(az_aspect)
az_matrix <- matrix(raster::extract(small_az,
raster::extent(small_az),buffer=500),
nrow=ncol(small_az), ncol=nrow(small_az))
az_matrix %>%
sphere_shade(texture = 'desert') %>%
add_water(detect_water(az_matrix),color='desert') %>%
add_shadow(ray_shade(az_matrix,zscale=2),0.5) %>%
add_shadow(ambient_shade(az_matrix)) %>%
plot_3d(az_matrix,zscale=10,fov = 0,theta=135,zoom=0.75)
render_snapshot(clear = TRUE)
For this section we are really going to focus on the mine area so first I’m going to buffer our az_mine polygon by 2km and then clip both the before and after mining rasters to this new shape
#Read in raster data
az_raster_after <- raster('data/srtm_14_06.tif')
az_mine_2km <- az_mine %>%
st_transform(2163) %>%
st_buffer(2000) %>%
#Reproject into SRTM datum (WGS84)
st_transform(st_crs(az_raster_after))
#Get extent of buffered mine area
aoi <- extent(az_mine_2km)
#Crop after to area of interest
mine_after <- crop(az_raster_after,aoi)
#Crop before to area of interest (need to also reproject)
mine_before <- crop(az_raster_before %>%
projectRaster(.,crs=projection(mine_after)),
aoi) %>%
projectRaster(.,mine_after) #match resolution
#Force rasters to be exactly same size
mine_before_same <- trim(mine_before)
mine_after_same <- crop(mine_after,mine_before_same)
#Stack these things in a rasterBrick (which we can do because we
#made them the exact same size)
mines <- brick(mine_before_same,mine_after_same)
#Add a dem_difference layer
mines[[3]] <- mines[[1]] - mines[[2]]
#Rename raster layers
names(mines) <- c('Pre_mining','Post_mining','DoD')
plot(mines)
All data for this section should be the data clipped to the mining area (generated in the code chunk above)
az_slope_after <- terrain(mine_after_same,opt='slope',unit='degrees')
mapview(az_slope_after)
az_slope_before <- terrain(mine_before_same,opt='slope',unit='degrees')
az_slope_diff<-terrain(mines[[3]], opt='slope', unit='degrees')
slope<-brick(az_slope_after,az_slope_before,az_slope_diff)
names(slope) <- c('Post_mining','Pre_mining','DoD')
plot(slope)
az_before_matrix <- matrix(raster::extract(mine_before_same,
raster::extent(mine_before_same)),
nrow=ncol(mine_before_same), ncol=nrow(mine_before_same))
az_before_matrix %>%
sphere_shade(texture = 'desert') %>%
add_water(detect_water(az_before_matrix),color='desert') %>%
add_shadow(ray_shade(az_before_matrix,zscale=2),0.5) %>%
add_shadow(ambient_shade(az_before_matrix)) %>%
plot_3d(az_before_matrix,zscale=10,fov = 0,theta=135,zoom=0.75)
render_snapshot(clear = TRUE)
az_after_matrix <- matrix(raster::extract(mine_after_same,
raster::extent(mine_after_same)),
nrow=ncol(mine_after_same), ncol=nrow(mine_after_same))
az_after_matrix %>%
sphere_shade(texture = 'desert') %>%
add_water(detect_water(az_after_matrix),color='desert') %>%
add_shadow(ray_shade(az_after_matrix,zscale=2),0.5) %>%
add_shadow(ambient_shade(az_after_matrix)) %>%
plot_3d(az_after_matrix,zscale=10,fov = 0,theta=135,zoom=0.75)
render_snapshot(clear = TRUE)